% ShapeFactorBinning
% Eden Ford
% 29 May 2023

% Code to calculate (1) the shape factor distribution and binning

% General definitions:
% Object - cell or cluster of cells that are touching/interacting with one
% another
% Shape factor - defines circularity (1 = perfect sphere)

close all
clear all
clc

%% 0 mM CMP1a Data Import
filename = 'Morphology.xlsx'; % file name to read
sheet = 1; % sheet to read

[num_0,txt_0,raw_0] = xlsread(filename,sheet);

ImageNames_0 = raw_0(2:end,2); % list of image each object is found in
ImageNumber_0 = unique(ImageNames_0); % create array of unique image names
Names_0 = char(ImageNames_0);
ImageNumberChar_0 = char(ImageNumber_0); % convert cell aray into character array
ShapeFactor_0 = cell2mat(raw_0(2:end,6)); % list of shape factor of each object

T_Full_0 = table(ImageNames_0,ShapeFactor_0) % organize all data into table

% Initialize shape factor counts
ShapeFactor_0_Gel1 = []; % cells per object count in 0 mM gel 1
ShapeFactor_0_Gel2 = []; % cells per object count in 0 mM gel 2
ShapeFactor_0_Gel3 = []; % cells per object count in 0 mM gel 3

for i = 1:length(Names_0)
    if Names_0(i,1:11) == '0mMCMP_Gel1'
        ShapeFactor_0_Gel1 = [ShapeFactor_0_Gel1,ShapeFactor_0(i)];
    elseif Names_0(i,1:11) == '0mMCMP_Gel2'
        ShapeFactor_0_Gel2 = [ShapeFactor_0_Gel2,ShapeFactor_0(i)];
    elseif Names_0(i,1:11) == '0mMCMP_Gel3'
        ShapeFactor_0_Gel3 = [ShapeFactor_0_Gel3,ShapeFactor_0(i)];
    end
end

ShapeFactor_0a = [ShapeFactor_0_Gel1 ShapeFactor_0_Gel2 ShapeFactor_0_Gel3]'; % should be equal to ShapeFactor_0

%% 10 mM CMP1a Data Import
filename = 'Morphology.xlsx'; % file name to read
sheet = 2; % sheet to read

[num_10,txt_10,raw_10] = xlsread(filename,sheet);

ImageNames_10 = raw_10(2:end,2); % list of image each object is found in
ImageNumber_10 = unique(ImageNames_10); % create array of unique image names
Names_10 = char(ImageNames_10);
ImageNumberChar_10 = char(ImageNumber_10); % convert cell aray into character array
ShapeFactor_10 = cell2mat(raw_10(2:end,6)); % list of shape factor of each object

T_Full_10 = table(ImageNames_10,ShapeFactor_10) % organize all data into table

% Initialize shape factor counts
ShapeFactor_10_Gel1 = []; % cells per object count in 10 mM gel 1
ShapeFactor_10_Gel2 = []; % cells per object count in 10 mM gel 2
ShapeFactor_10_Gel3 = []; % cells per object count in 10 mM gel 3

for i = 1:length(Names_10)
    if Names_10(i,1:12) == '10mMCMP_Gel1'
        ShapeFactor_10_Gel1 = [ShapeFactor_10_Gel1,ShapeFactor_10(i)];
    elseif Names_10(i,1:12) == '10mMCMP_Gel2'
        ShapeFactor_10_Gel2 = [ShapeFactor_10_Gel2,ShapeFactor_10(i)];
    elseif Names_10(i,1:12) == '10mMCMP_Gel3'
        ShapeFactor_10_Gel3 = [ShapeFactor_10_Gel3,ShapeFactor_10(i)];
    end
end

ShapeFactor_10a = [ShapeFactor_10_Gel1 ShapeFactor_10_Gel2 ShapeFactor_10_Gel3]'; % should be equal to ShapeFactor_0

%% Shape Factor Distribution and Statistical Analysis

%ShapeFactor_0
%ShapeFactor_10

% Display results
disp(' ');
disp('Shape factor analysis');
disp(' ');

% Moments of distribution
Type_ShapeFactor = {'0mMCMP'; '10mMCMP'};
% MEAN
Average_ShapeFactor = [mean(ShapeFactor_0); mean(ShapeFactor_10)];
% STANDARD ERROR (DERIVED FROM VARIANCE)
StError_ShapeFactor = [std(ShapeFactor_0/sqrt(3)); std(ShapeFactor_10/sqrt(3))];
% MEDIAN
Median_ShapeFactor = [median(ShapeFactor_0); median(ShapeFactor_10)];
% MODE
Mode_ShapeFactor = [mode(ShapeFactor_0); mode(ShapeFactor_10)];
% SKEWNESS
Skewness_ShapeFactor = [skewness(ShapeFactor_0); skewness(ShapeFactor_10)];
% KURTOSIS
Kurtosis_ShapeFactor = [kurtosis(ShapeFactor_0); kurtosis(ShapeFactor_10)];

ShapeFactorTable = table(Type_ShapeFactor, Average_ShapeFactor, StError_ShapeFactor, Median_ShapeFactor, Mode_ShapeFactor, Skewness_ShapeFactor, Kurtosis_ShapeFactor);
ShapeFactorTable.Properties.VariableNames = {'GelComposition' 'Average' 'StError' 'Median' 'Mode' 'Skewness' 'Kurtosis'};
disp(ShapeFactorTable);

% Statistical analysis
filename = 'Morphology.xlsx'; % file name to read
sheet = 3; % sheet to read

[num,txt,raw] = xlsread(filename,sheet);
ShapeFactorFull = num(3:end,:);

[p tbl stats] = anova1(ShapeFactorFull) % where T is a matrix with data columns
ShapeFactor_comparison = multcompare(stats) % multiple comparisons test[h2a_RL3,p2a_RL3] = ttest2(ShapeFactor_0mMCMP,ShapeFactor2a_Low);

%% Histograms

figure
BinWidth = 0.05;
yMax = 0.3;
edges = linspace(0, 1, 21);
histogram_distribution = zeros(6,20);

    % 0 mM CMP
    H2a = subplot(2,1,1);
    hold on
    histogram(ShapeFactor_0, 'Normalization', 'probability', 'BinWidth', BinWidth);
    ylabel('Fraction of Objects')
    hold off
    title('0 mM CMP')
    xlim([0 1])
    ylim([0 yMax])
    histogram_distribution(1,:) = histcounts(ShapeFactor_0,edges);
    
    % 10 mM CMP
    H2b = subplot(2,1,2);
    hold on
    histogram(ShapeFactor_10, 'Normalization', 'probability', 'BinWidth', BinWidth);
    % ylabel('Fraction of Objects')
    hold off
    title('10 mM CMP')
    xlim([0 1])
    ylim([0 yMax])
    histogram_distribution(2,:) = histcounts(ShapeFactor_10,edges);

%% Binning

% Bin1 = Low (>0 - =0.275)
% Bin2 = Middle (>0.275 - =0.475)
% Bin3 = High (>0.475 - =1)
Cutoff1 = 0.25;
Cutoff2 = 0.5;
BinName = ['gel1,rounded ','gel1,spread ','gel1,elongated';'gel2,rounded ','gel2,spread ','gel2,elongated';'gel3,rounded ','gel3,spread ','gel3,elongated']

% 0 mM CMP
% Initialize bins
Bin = [0 0 0; 0 0 0; 0 0 0];
for i = 1:length(ShapeFactor_0_Gel1) % 0 mM CMP, Gel 1
    if ShapeFactor_0_Gel1(i) > Cutoff2
        Bin(1,1) = Bin(1,1) + 1;
    elseif ShapeFactor_0_Gel1(i) > Cutoff1
        Bin(1,2) = Bin(1,2) +1;
    elseif ShapeFactor_0_Gel1(i) > 0
        Bin(1,3) = Bin(1,3) +1;
    end
end

                                                                                                                                                        
for i = 1:length(ShapeFactor_0_Gel2) % 0 mM CMP, Gel 2
    if ShapeFactor_0_Gel2(i) > Cutoff2
        Bin(2,1) = Bin(2,1) + 1;
    elseif ShapeFactor_0_Gel2(i) > Cutoff1
        Bin(2,2) = Bin(2,2) +1;
    elseif ShapeFactor_0_Gel2(i) > 0
        Bin(2,3) = Bin(2,3) +1;
    end
end

for i = 1:length(ShapeFactor_0_Gel3) % 0 mM CMP, Gel 3
    if ShapeFactor_0_Gel3(i) > Cutoff2
        Bin(3,1) = Bin(3,1) + 1;
    elseif ShapeFactor_0_Gel3(i) > Cutoff1
        Bin(3,2) = Bin(3,2) +1;
    elseif ShapeFactor_0_Gel3(i) > 0                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                
        Bin(3,3) = Bin(3,3) +1;                                                                                 
    end
end

ShapeFactor_0_Count = Bin

% 10 mM CMP
% Initialize bins
Bin = [0 0 0; 0 0 0; 0 0 0];
for i = 1:length(ShapeFactor_10_Gel1) % 10 mM CMP, Gel 1
    if ShapeFactor_10_Gel1(i) > Cutoff2
        Bin(1,1) = Bin(1,1) + 1;
    elseif ShapeFactor_10_Gel1(i) > Cutoff1
        Bin(1,2) = Bin(1,2) +1;
    elseif ShapeFactor_10_Gel1(i) > 0
        Bin(1,3) = Bin(1,3) +1;
    end
end
                                                                                                                                                        
for i = 1:length(ShapeFactor_10_Gel2) % 10 mM CMP, Gel 2
    if ShapeFactor_10_Gel2(i) > Cutoff2
        Bin(2,1) = Bin(2,1) + 1;
    elseif ShapeFactor_10_Gel2(i) > Cutoff1
        Bin(2,2) = Bin(2,2) +1;
    elseif ShapeFactor_10_Gel2(i) > 0
        Bin(2,3) = Bin(2,3) +1;
    end
end

for i = 1:length(ShapeFactor_10_Gel3) % 10 mM CMP, Gel 3
    if ShapeFactor_10_Gel3(i) > Cutoff2
        Bin(3,1) = Bin(3,1) + 1;
    elseif ShapeFactor_10_Gel3(i) > Cutoff1
        Bin(3,2) = Bin(3,2) +1;
    elseif ShapeFactor_10_Gel3(i) > 0                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                
        Bin(3,3) = Bin(3,3) +1;                                                                                 
    end
end

ShapeFactor_10_Count = Bin